# SLE human patients
source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

# examine mva pathway ----------
kegg_mva <-
  c('Acat1','Acat2','Hmgcs1','Hmgcs2','Hmgcr','Mvk','Pmvk','Mvd','Idi1','Idi2','Fdps') |>
  str_to_upper()

mva_fct <- tibble(gene = kegg_mva, ordered = fct_inorder(kegg_mva))

# examine interested cytokines ------
key_cytokine <- c('Il6','Ccl20','Tslp','Flt3lg','Csf2','Tnf') |>
  str_to_upper()

# read mtx data ----------
safe10x <- safely(Read10X, quiet = F)

# turn warning to error
options(warn = 2)
## zheng2022 --------
zh22_dir <- list.dirs('mission/FPP/xiangya_sle_scRNA/data/zheng2022/',
                      full.names = T, recursive = F)

zh22_dir <- zh22_dir |>
  basename() |>
  set_names(zh22_dir, nm = _)

# 2 batches have different genes list, need use intersection
batch1_zh22 <- zh22_dir[1:8] |>
  Read10X(strip.suffix = T)

batch2_zh22 <- zh22_dir[9:17] |>
  Read10X(strip.suffix = T)

zh22b1_name <- zh22_dir[1:8] |> names()

batch1_zh22[1:5,1:5]

batch2_zh22[1:5,1:5]

batch1_zh22 |> glimpse()
batch2_zh22 |> glimpse()

intsc_genes <- rownames(batch1_zh22) |>
  intersect(rownames(batch2_zh22))

## new xiangya data -----
nxy_dir <- list.dirs('mission/FPP/xiangya/',
                      full.names = T, recursive = F) |>
  str_subset('SLE')

nxy_dir <- nxy_dir |>
  basename() |>
  str_c('New') |>
  set_names(nxy_dir, nm = _)

batch_nxy <- nxy_dir |>
  Read10X(strip.suffix = T)

batch_nxy |> glimpse()

intsc_genes <- rownames(batch_nxy) |>
  intersect(intsc_genes)

## kim2023 --------
kim_dir <- list.dirs('mission/FPP/xiangya_sle_scRNA/data/kim2023/',
                     full.names = T, recursive = F)

kim_dir <- kim_dir |>
  basename() |>
  str_c('kim') |>
  set_names(kim_dir, nm = _)

batch_kim <- kim_dir |>
  Read10X(strip.suffix = T)

batch_kim |> glimpse()

intsc_genes <- rownames(batch_kim) |>
  intersect(intsc_genes)

# turn to normal warn level
options(warn = 1)

# create sobj -------
sobj <- list(batch1_zh22, batch2_zh22) |>
  map(\(x)x[intsc_genes,], .progress = T) |>
  purrr::reduce(RowMergeSparseMatrices) |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

sobj$mito.ratio <- sobj |>
  PercentageFeatureSet('^MT-')

sobj |> VlnPlot('mito.ratio', pt.size = 0)

sobj <- sobj |>
  filter(mito.ratio < 10)

sobj <- sobj |>
  mutate(sle = str_extract(orig.ident, 'HC|SLE'),
         dataset = case_when(
           str_detect(orig.ident, 'kim') ~ 'kim2023',
           str_detect(orig.ident, 'New') ~ 'xiangya_new',
           orig.ident %in% zh22b1_name ~ 'zheng2022_batch1',
           .default = 'zheng2022_batch2'
         ),
         tissue = ifelse(str_detect(orig.ident, 'D$'), 'dermis', 'epidermis'))

sobj |>
  VlnPlot('mito.ratio', group.by = 'tissue', pt.size = 0)

sobj <- sobj |>
  quick_process_seurat(batch = c('orig.ident', 'dataset', 'tissue'))

sobj |> DimPlot(group.by = 'dataset')

sobj |> DimPlot(group.by = 'sle')

## auto annotate cell type ---------
hpca <- celldex::HumanPrimaryCellAtlasData()

sobj <- sobj |>
  mark_cell_type_singler(ref = hpca, new_label = 'hpca_main')

sobj |> DimPlot(group.by = 'hpca_main', cols = DiscretePalette(36),
                label = T, label.box = T)

sobj |>
  ggplot(aes(tissue, fill = hpca_main)) +
  geom_bar(position = 'fill') +
  scale_fill_manual(values = DiscretePalette(36))

sobj |>
  ggplot(aes(sle, fill = hpca_main)) +
  geom_bar(position = 'fill') +
  scale_fill_manual(values = DiscretePalette(36))

sobj |>
  ggplot(aes(tissue, fill = hpca_main)) +
  geom_bar(position = 'fill') +
  scale_fill_manual(values = DiscretePalette(36)) +
  theme_pubr()

sobj |>
  dplyr::count(hpca_main, sle, tissue) |>
  DT::datatable()

## TRPM2 expr -------
sobj |>
  DotPlot2d('TRPM2', sle, hpca_main)

sobj |>
  mutate(sle_tissue = str_c(sle, '_', tissue)) |>
  DotPlot2d('TRPM2', sle_tissue, hpca_main) +
  labs(title = 'TRPM2 expression in SLE skin tissue',
       x = 'Sample', y = 'Cell type')

g1 <- last_plot()

g1 + RotatedAxis()

# keratinocyte & melanocyte marker genes
sobj |> DotPlot(c('KRT1','KRT14','PMEL','MLANA'), cluster.idents = T) +
  coord_flip()

sobj |>
  dplyr::count(orig.ident, individual, dataset) |>
  write_tsv('human_trpv3_sle.sample.tsv')

## manual edit cell type -----
sobj <- sobj |>
  mutate(manual_main = case_when(seurat_clusters == 22 ~ 'Melanocytes',
                                 hpca_main == 'Tissue_stem_cells' ~ 'Fibroblasts',
                                 hpca_main == 'Epithelial_cells' ~ 'Keratinocytes',
                                 .default = hpca_main))

sobj |> DimPlot(group.by = 'manual_main', cols = DiscretePalette(36),
                label = T, label.box = T, repel = T) +
  ggtitle('Cell types in SLE dermis + epidermis')

# select only epidermis ----
sobj_epi <- sobj |>
  filter(tissue == 'epidermis')

sobj_epi |> write_rds('mission/FPP/human_SLE_epi.rds')

# find trpv3h in kim2023 -----------
sobj_kkera <- sobj |>
  filter(dataset == 'kim2023' & hpca_main == 'Keratinocytes')

sobj_kkera <- sobj_kkera |>
  quick_process_seurat(skip_norm = T)

sobj_kkera |> DotPlot('TRPV3')

sobj_kkera |> FindAllMarkers(features = 'TRPV3',
                             only.pos = T)

sobj_kkera |> summarise(n(), .by = seurat_clusters)

# CCR6 in cDC2? --------
sobj.epi <- read_rds('mission/FPP/xiangya_sle_scRNA/data/human_SLE_epi.rds')

sle.dc <- sobj.epi |> filter(bill_fine == 'DC')
sle.dc |> FeaturePlot('CCR6')

sle.dc$orig.ident |> table()

sle.dc <- sle.dc |>
  quick_process_seurat(skip_norm = T)

sle.dc |>
  DotPlot(c('CD1C','FCER1A',
            'SIRPA','CCR6'), cluster.idents = T)

sle.dc <- sle.dc |>
  mutate(dc.subtype = ifelse(seurat_clusters %in% c(11,5,7), 'cDC2',
                             'other DC')) 

sle.dc |>
  DimPlot(group.by = 'dc.subtype')

sle.dc |>
  VlnPlot('CCR6', group.by = 'dc.subtype')

sle.dc |>
  FindMarkers(features = 'CCR6', group.by = 'dc.subtype', ident.1 = 'cDC2')

## TRPM2-hi macrophage -----------
sobj.epi |>
  DotPlot2d('TRPM2', sle, bill_fine) +
  labs(x = 'Group', y = 'Cell type',
       title = 'TRPM2 expression in SLE epidermis')

# check new 3v4 mixed skin -------
new.path <- list.files('mission/FPP/xiangya_sle_scRNA/data',
                       include.dirs = T, pattern = 'SLE\\d',
                       full.names = T)

nsle.mex <- new.path |>
  str_extract('SLE\\d') |>
  set_names(x = new.path, nm = _) |>
  Read10X()

nsle.mex |> glimpse()

new.hc <- list.files('mission/FPP/xiangya_sle_scRNA/data',
                       include.dirs = T, pattern = 'HC\\d',
                       full.names = T)

nhc.mex <- new.hc |>
  str_extract('HC\\d') |>
  set_names(x = new.hc, nm = _) |>
  Read10X()

nhc.mex |> glimpse()

nset.gene <- rownames(nsle.mex) |>
  intersect(rownames(nhc.mex))

length(nset.gene)

sobj <- nsle.mex[nset.gene, ] |>
  RowMergeSparseMatrices(nhc.mex[nset.gene, ]) |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

sobj$mito.ratio <- sobj |>
  PercentageFeatureSet('^MT-')

sobj |>
  VlnPlot('mito.ratio', pt.size = 0)

sobj <- sobj |>
  filter(mito.ratio < 10)

sobj <- sobj |>
  mutate(group = str_extract(orig.ident, 'HC|SLE'))

sobj |> VlnPlot(c('nCount_RNA', 'nFeature_RNA'),
                group.by = 'group', pt.size = 0)

sobj |>
  ggplot(aes(group, log10(nFeature_RNA))) +
  geom_boxplot()

sobj <- sobj |>
  quick_process_seurat(batch = c('group','orig.ident'))

hpca <- celldex::HumanPrimaryCellAtlasData()

sobj <- sobj |>
  mark_cell_type_singler(ref = hpca, new_label = 'hpca.main')

sobj |>
  DimPlot(group.by = 'hpca.main', cols = 'Paired')

sobj |>
  DotPlot(c('KRT1','KRT14'), cluster.idents = T,
          cols = 'RdYlBu')

sobj |>
  write_rds('mission/FPP/xiangya_sle_scRNA/new.mixed.rds')

sobj <- read_rds('mission/FPP/xiangya_sle_scRNA/new.mixed.rds')

kera.am <- sobj |>
  FindAllMarkers(features = c('KRT1','KRT14'), only.pos = T) |>
  filter(p_val_adj < .05) |>
  pull(cluster)

sobj |>
  DimPlot(cols = DiscretePalette(36), label = T, label.box = T)

sobj <- sobj |>
  mutate(manual.main = case_when(seurat_clusters %in% c(21,15,14,16,17,19) ~ 'Keratinocytes',
                                 .default = hpca.main))

### UMAP cell type ----------
sobj |>
  DimPlot(group.by = 'manual.main', cols = DiscretePalette(36),
          label = T, label.size = 2, repel = T) +
  ggtitle('Cell types of skin') +
  theme_jpub +
  NoLegend()

publish_pdf('mission/FPP/figures/human.mixed.sle.celltype.umap.pdf')

sobj |>
  filter(group == 'HC') |>
  DotPlot(key_cytokine, cols = 'RdYlBu', group.by = 'manual.main')

apf.main.expr <- last_plot() |>
  pluck('data')

apf.main.expr |>
  ggplot(aes(features.plot, id, fill = avg.exp.scaled, size = pct.exp)) +
  geom_point(shape = 21, stroke = .3) +
  scale_fill_gradient2(low = 'blue', high = 'red') +
  theme_bw() +
  #theme_jpub +
  RotatedAxis() +
  labs(title = 'Human HC mixed epidermis & dermis', x = 'Gene', y = 'Cell type',
       fill = 'Average expression', size = 'Percent expressed')

## KC ---------
sobj.kcn <- sobj |>
  filter(str_detect(manual.main, 'Kera'))

sobj.kcn <- sobj.kcn |>
  quick_process_seurat(batch = c('orig.ident', 'group'))

sobj.kcn <- sobj.kcn |>
  mark_cell_type_singler(ref = hpca, new_label = 'hpca.main2')

sobj.kcn |>
  DimPlot(group.by = 'hpca.main2')

### V3 high/low -----------
late.kc <- c('Krt1','Krt10','Lor','Ivl','Tgm1','Flg') |>
  str_to_upper()

g1 <- sobj.kcn |>
  DotPlot(c('TRPV3', late.kc, 'KRT5', 'KRT14'), cluster.idents = T)

mixed.kc.mark <- g1 |>
  pluck('data')

mixed.kc.mark |>
  mutate(id = as.character(id) |> fct_relevel('2','6','7','11')) |>
  ggplot(aes(features.plot, id, size = pct.exp)) +
  scale_fill_gradient2(low = 'blue', high = 'red') +
  theme_pubr(legend = "right") +
  labs(fill = "Average expression", size = "Percent expressed",
       x = 'Gene', y = 'Cell type') +
  geom_point(aes(fill = avg.exp.scaled), shape = 21, stroke = .3) +
  theme_jpub +
  RotatedAxis() +
  scale_size(range = c(0,2))

publish_pdf('mission/FPP/figures/mixed.kc.v3.early.late.marker.pdf',
            width = 70)

sobj.kcn |>
  FindAllMarkers(features = 'TRPV3', only.pos = T)

sobj.kcn <- sobj.kcn |>
  mutate(v3.fine = case_when(seurat_clusters == 21 ~ 'T_cells',
                             seurat_clusters %in% c(2,6,7,11) ~ 'Trpv3-hi-KC',
                             .default = 'Trpv3-lo-KC'))

sobj.kcn |>
  DotPlot2d('CCL20', v3.fine, group)

sobj <- sobj.kcn |>
  as_tibble() |>
  select(.cell, v3.fine) |>
  left_join(x = sobj, y = _)
  
sobj <- sobj |>
  mutate(manual.fine = case_when(is.na(v3.fine) ~ manual.main,
                                 .default = v3.fine))
### MVA expr in SLE -------------
sobj |>
  filter(group == 'SLE') |>
  BubblePlot(kegg_mva, group.by = 'manual.fine') +
  RotatedAxis()

g1 <- last_plot()

sle.mva.expr <- g1 |>
  pluck('data')

g1 + scale_size(range = c(0, 4)) +
  theme_jpub +
  RotatedAxis()

publish_pdf('mission/FPP/figures/mixed.sle.skin.mva.expr.pdf', width = 70)

sobj |>
  filter(group == 'SLE') |>
  BubblePlot(kegg_mva, group.by = 'manual.fine') +
  RotatedAxis()

### APF expr in SLE -------------
sobj |>
  filter(group == 'SLE') |>
  BubblePlot(key_cytokine, group.by = 'manual.fine') +
  RotatedAxis()

g1 <- last_plot()

sle.apf.expr <- g1 |>
  pluck('data')

g1 + scale_size(range = c(0, 5)) +
  theme_jpub +
  RotatedAxis() + labs(x = 'Gene', y = 'Cell type')

publish_pdf('mission/FPP/figures/mixed.sle.skin.apf.expr.pdf', width = 60)

### KC diff marker expr -------------
sobj |>
  BubblePlot(c('TRPV3', late.kc, 'KRT5', 'KRT14'),
             group.by = 'manual.fine') +
  RotatedAxis()

g1 <- last_plot()

mixed.kcdiff.expr <- g1 |>
  pluck('data')

g1 + scale_size(range = c(0,3)) +
  theme_jpub +
  RotatedAxis() + labs(x = 'Gene', y = 'Cell type')

publish_pdf('mission/FPP/figures/mixed.skin.kcdiff.expr.pdf', width = 65)

### SLE vs HC DEGs ----------
mixed.slevhc.deg <- sobj |>
  FindMarkersAcrossVar(split.by = 'manual.fine', group.by = 'group',
                       ident.1 = 'SLE', logfc.threshold = .01, min.pct = 0)

mixed.slevhc.deg |>
  write_csv('mission/FPP/xiangya_sle_scRNA/results/mixed.slevhc.deg.csv')

mixed.slevhc.deg |>
  filter(cluster == 'Trpv3-hi-KC') |>
  plot_bill_volc(group1 = 'SLE v3h', group2 = 'HC v3h',
                 highlights = c('HMGCR', 'HMCGS1', 'MVK', key_cytokine),
                 force = T)

mixed.slevhc.deg |>
  right_join(mva_fct) |>
  ggplot(aes(ordered, cluster, fill = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point(shape = 21, stroke = .3) +
  scale_fill_gradient2(low = 'blue', high = 'red') +
  scale_size(range = c(0,3)) +
  theme_jpub +
  RotatedAxis() +
  labs(x = 'Gene', y = 'Cell type')

publish_pdf('mission/FPP/figures/mixed.slevhc.logfc.mva.pdf', width = 70)

mixed.slevhc.deg |> 
  filter(gene %in% key_cytokine) |>
  mutate(p_val_adj = ifelse(p_val_adj == 0, 1e-300, p_val_adj)) |>
  ggplot(aes(gene, cluster, fill = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point(shape = 21, stroke = .3) +
  scale_fill_gradient2(low = 'blue', high = 'red') +
  scale_size(range = c(0,3)) +
  theme_jpub +
  RotatedAxis() +
  ylab('Cell type')
 
publish_pdf('mission/FPP/figures/mixed.slevhc.logfc.apf.pdf', width = 60)

mixed.slevhc.deg.kc <- sobj |>
  filter(str_detect(manual.fine, 'KC')) |>
  FindMarkers(group.by = 'group', ident.1 = 'SLE',
              logfc.threshold = .01, min.pct = 0) |>
  as_tibble(rownames = 'gene')

mixed.slevhc.deg.kc |>
  plot_bill_volc(group1 = 'SLE KC', group2 = 'HC KC',
                 highlights = c('HMGCR', 'HMCGS1', 'MVK', key_cytokine),
                 force = T)

### V3-hi vs V3-lo in SLE ----------
sle.v3hvl.deg <- sobj.kcn |>
  filter(group == 'SLE') |>
  FindMarkers(group.by = 'v3.fine',
              ident.1 = 'Trpv3-hi-KC', ident.2 = 'Trpv3-lo-KC')

sle.v3hvl.deg |>
  plot_pub_volc(group1 = 'SLE TRPV3-hi-KC', group2 = 'SLE TRPV3-lo-KC',
                 highlights = c('HMGCR', 'HMCGS1', 'MVK', key_cytokine),
                 force = T)

sle.v3hvl.deg |>
  write_csv('mission/FPP/xiangya_sle_scRNA/results/mixed.sle.v3hvl.deg.csv')

## tpm percent ------------
ccl20.tpm.hm <- sobj |>
  get_abundance_sc_wide('CCL20') |>
  mutate(CCL20 = expm1(CCL20)) |>
  left_join(x = sobj, y = _) |>
  summarise(sum = sum(CCL20), .by = c(manual.main, group))

ccl20.tpm.hm |>
  filter(group == 'HC') |>
  mutate(percent = sum / sum(sum))

# ADRB2 / ADRA1A ----------
sobj.epi |>
  DotPlot(c('ADRB2','ADRA1A'))
